1//////////////////////////////////////////////////////////////////////
   2// LibFile: paths.scad
   3//   A `path` is a list of points of the same dimensions, usually 2D or 3D, that can
   4//   be connected together to form a sequence of line segments or a polygon.
   5//   A `region` is a list of paths that represent polygons, and the functions
   6//   in this file work on paths and also 1-regions, which are regions
   7//   that include exactly one path.  When you pass a 1-region to a function, the default
   8//   value for `closed` is always `true` because regions represent polygons.  
   9//   Capabilities include computing length of paths, computing
  10//   path tangents and normals, resampling of paths, and cutting paths up into smaller paths.  
  11// Includes:
  12//   include <BOSL2/std.scad>
  13// FileGroup: Advanced Modeling
  14// FileSummary: Operations on paths: length, resampling, tangents, splitting into subpaths
  15// FileFootnotes: STD=Included in std.scad
  16//////////////////////////////////////////////////////////////////////
  17
  18// Section: Utility Functions
  19
  20// Function: is_path()
  21// Synopsis: Returns True if 'list' is a path.
  22// Topics: Paths
  23// See Also: is_region(), is_vnf()
  24// Usage:
  25//   is_path(list, [dim], [fast])
  26// Description:
  27//   Returns true if `list` is a path.  A path is a list of two or more numeric vectors (AKA points).
  28//   All vectors must of the same size, and may only contain numbers that are not inf or nan.
  29//   By default the vectors in a path must be 2d or 3d.  Set the `dim` parameter to specify a list
  30//   of allowed dimensions, or set it to `undef` to allow any dimension.  (Note that this function
  31//   returns `false` on 1-regions.)  
  32// Example:
  33//   bool1 = is_path([[3,4],[5,6]]);    // Returns true
  34//   bool2 = is_path([[3,4]]);          // Returns false
  35//   bool3 = is_path([[3,4],[4,5]],2);  // Returns true
  36//   bool4 = is_path([[3,4,3],[5,4,5]],2);  // Returns false
  37//   bool5 = is_path([[3,4,3],[5,4,5]],2);  // Returns false
  38//   bool6 = is_path([[3,4,5],undef,[4,5,6]]);  // Returns false
  39//   bool7 = is_path([[3,5],[undef,undef],[4,5]]);  // Returns false
  40//   bool8 = is_path([[3,4],[5,6],[5,3]]);     // Returns true
  41//   bool9 = is_path([3,4,5,6,7,8]);           // Returns false
  42//   bool10 = is_path([[3,4],[5,6]], dim=[2,3]);// Returns true
  43//   bool11 = is_path([[3,4],[5,6]], dim=[1,3]);// Returns false
  44//   bool12 = is_path([[3,4],"hello"], fast=true); // Returns true
  45//   bool13 = is_path([[3,4],[3,4,5]]);            // Returns false
  46//   bool14 = is_path([[1,2,3,4],[2,3,4,5]]);      // Returns false
  47//   bool15 = is_path([[1,2,3,4],[2,3,4,5]],undef);// Returns true
  48// Arguments:
  49//   list = list to check
  50//   dim = list of allowed dimensions of the vectors in the path.  Default: [2,3]
  51//   fast = set to true for fast check that only looks at first entry.  Default: false
  52function is_path(list, dim=[2,3], fast=false) =
  53    fast
  54    ?   is_list(list) && is_vector(list[0]) 
  55    :   is_matrix(list) 
  56        && len(list)>1 
  57        && len(list[0])>0
  58        && (is_undef(dim) || in_list(len(list[0]), force_list(dim)));
  59
  60// Function: is_1region()
  61// Synopsis: Returns true if path is a region with one component.
  62// Topics: Paths, Regions
  63// See Also: force_path()
  64// Usage:
  65//   bool = is_1region(path, [name])
  66// Description:
  67//   If `path` is a region with one component (a 1-region) then return true.  If path is a region with more components
  68//   then display an error message about the parameter `name` requiring a path or a single component region.  If the input
  69//   is not a region then return false.  This function helps path functions accept 1-regions.
  70// Arguments:
  71//   path = input to process
  72//   name = name of parameter to use in error message.  Default: "path"
  73function is_1region(path, name="path") = 
  74     !is_region(path)? false
  75    :assert(len(path)==1,str("Parameter \"",name,"\" must be a path or singleton region, but is a multicomponent region"))
  76     true;
  77
  78
  79// Function: force_path()
  80// Synopsis: Checks that path is a region with one component.
  81// SynTags: Path
  82// Topics: Paths, Regions
  83// See Also: is_1region()
  84// Usage:
  85//   outpath = force_path(path, [name])
  86// Description:
  87//   If `path` is a region with one component (a 1-region) then returns that component as a path.
  88//   If path is a region with more components then displays an error message about the parameter
  89//   `name` requiring a path or a single component region.  If the input is not a region then
  90//   returns the input without any checks.  This function helps path functions accept 1-regions.
  91// Arguments:
  92//   path = input to process
  93//   name = name of parameter to use in error message.  Default: "path"
  94function force_path(path, name="path") =
  95   is_region(path) ?
  96       assert(len(path)==1, str("Parameter \"",name,"\" must be a path or singleton region, but is a multicomponent region"))
  97       path[0]
  98   : path;
  99
 100
 101/// Internal Function: _path_select()
 102/// Usage:
 103///   _path_select(path,s1,u1,s2,u2,[closed]):
 104/// Description:
 105///   Returns a portion of a path, from between the `u1` part of segment `s1`, to the `u2` part of
 106///   segment `s2`.  Both `u1` and `u2` are values between 0.0 and 1.0, inclusive, where 0 is the start
 107///   of the segment, and 1 is the end.  Both `s1` and `s2` are integers, where 0 is the first segment.
 108/// Arguments:
 109///   path = The path to get a section of.
 110///   s1 = The number of the starting segment.
 111///   u1 = The proportion along the starting segment, between 0.0 and 1.0, inclusive.
 112///   s2 = The number of the ending segment.
 113///   u2 = The proportion along the ending segment, between 0.0 and 1.0, inclusive.
 114///   closed = If true, treat path as a closed polygon.
 115function _path_select(path, s1, u1, s2, u2, closed=false) =
 116    let(
 117        lp = len(path),
 118        l = lp-(closed?0:1),
 119        u1 = s1<0? 0 : s1>l? 1 : u1,
 120        u2 = s2<0? 0 : s2>l? 1 : u2,
 121        s1 = constrain(s1,0,l),
 122        s2 = constrain(s2,0,l),
 123        pathout = concat(
 124            (s1<l && u1<1)? [lerp(path[s1],path[(s1+1)%lp],u1)] : [],
 125            [for (i=[s1+1:1:s2]) path[i]],
 126            (s2<l && u2>0)? [lerp(path[s2],path[(s2+1)%lp],u2)] : []
 127        )
 128    ) pathout;
 129
 130
 131
 132// Function: path_merge_collinear()
 133// Synopsis: Removes unnecessary points from a path.
 134// SynTags: Path
 135// Topics: Paths, Regions
 136// Description:
 137//   Takes a path and removes unnecessary sequential collinear points.  Note that when `closed=true` either of the path
 138//   endpoints may be removed.  
 139// Usage:
 140//   path_merge_collinear(path, [eps])
 141// Arguments:
 142//   path = A path of any dimension or a 1-region
 143//   closed = treat as closed polygon.  Default: false
 144//   eps = Largest positional variance allowed.  Default: `EPSILON` (1-e9)
 145function path_merge_collinear(path, closed, eps=EPSILON) =
 146    is_1region(path) ? path_merge_collinear(path[0], default(closed,true), eps) :
 147    let(closed=default(closed,false))
 148    assert(is_bool(closed))
 149    assert( is_path(path), "Invalid path in path_merge_collinear." )
 150    assert( is_undef(eps) || (is_finite(eps) && (eps>=0) ), "Invalid tolerance." )
 151    len(path)<=2 ? path :
 152    let(path = deduplicate(path, closed=closed))
 153    [
 154      if(!closed) path[0],
 155      for(triple=triplet(path,wrap=closed))
 156        if (!is_collinear(triple,eps=eps)) triple[1],
 157      if(!closed) last(path)
 158    ];
 159
 160
 161// Section: Path length calculation
 162
 163
 164// Function: path_length()
 165// Synopsis: Returns the path length.
 166// Topics: Paths
 167// See Also: path_segment_lengths(), path_length_fractions()
 168// Usage:
 169//   path_length(path,[closed])
 170// Description:
 171//   Returns the length of the path.
 172// Arguments:
 173//   path = Path of any dimension or 1-region. 
 174//   closed = true if the path is closed.  Default: false
 175// Example:
 176//   path = [[0,0], [5,35], [60,-25], [80,0]];
 177//   echo(path_length(path));
 178function path_length(path,closed) =
 179    is_1region(path) ? path_length(path[0], default(closed,true)) :
 180    assert(is_path(path), "Invalid path in path_length")
 181    let(closed=default(closed,false))
 182    assert(is_bool(closed))
 183    len(path)<2? 0 :
 184    sum([for (i = [0:1:len(path)-2]) norm(path[i+1]-path[i])])+(closed?norm(path[len(path)-1]-path[0]):0);
 185
 186
 187// Function: path_segment_lengths()
 188// Synopsis: Returns a list of the lengths of segments in a path.
 189// Topics: Paths
 190// See Also: path_length(), path_length_fractions()
 191// Usage:
 192//   path_segment_lengths(path,[closed])
 193// Description:
 194//   Returns list of the length of each segment in a path
 195// Arguments:
 196//   path = path in any dimension or 1-region
 197//   closed = true if the path is closed.  Default: false
 198function path_segment_lengths(path, closed) =
 199    is_1region(path) ? path_segment_lengths(path[0], default(closed,true)) :
 200    let(closed=default(closed,false))
 201    assert(is_path(path),"Invalid path in path_segment_lengths.")
 202    assert(is_bool(closed))
 203    [
 204        for (i=[0:1:len(path)-2]) norm(path[i+1]-path[i]),
 205        if (closed) norm(path[0]-last(path))
 206    ]; 
 207
 208
 209// Function: path_length_fractions()
 210// Synopsis: Returns the fractional distance of each point along the length of a path.
 211// Topics: Paths
 212// See Also: path_length(), path_segment_lengths()
 213// Usage:
 214//   fracs = path_length_fractions(path, [closed]);
 215// Description:
 216//    Returns the distance fraction of each point in the path along the path, so the first
 217//    point is zero and the final point is 1.  If the path is closed the length of the output
 218//    will have one extra point because of the final connecting segment that connects the last
 219//    point of the path to the first point.
 220// Arguments:
 221//    path = path in any dimension or a 1-region
 222//    closed = set to true if path is closed.  Default: false
 223function path_length_fractions(path, closed) =
 224    is_1region(path) ? path_length_fractions(path[0], default(closed,true)):
 225    let(closed=default(closed, false))
 226    assert(is_path(path))
 227    assert(is_bool(closed))
 228    let(
 229        lengths = [
 230            0,
 231            each path_segment_lengths(path,closed)
 232        ],
 233        partial_len = cumsum(lengths),
 234        total_len = last(partial_len)
 235    )
 236    partial_len / total_len;
 237
 238
 239
 240/// Internal Function: _path_self_intersections()
 241/// Usage:
 242///   isects = _path_self_intersections(path, [closed], [eps]);
 243/// Description:
 244///   Locates all self intersection points of the given path.  Returns a list of intersections, where
 245///   each intersection is a list like [POINT, SEGNUM1, PROPORTION1, SEGNUM2, PROPORTION2] where
 246///   POINT is the coordinates of the intersection point, SEGNUMs are the integer indices of the
 247///   intersecting segments along the path, and the PROPORTIONS are the 0.0 to 1.0 proportions
 248///   of how far along those segments they intersect at.  A proportion of 0.0 indicates the start
 249///   of the segment, and a proportion of 1.0 indicates the end of the segment.
 250///   .
 251///   Note that this function does not return self-intersecting segments, only the points
 252///   where non-parallel segments intersect.  
 253/// Arguments:
 254///   path = The path to find self intersections of.
 255///   closed = If true, treat path like a closed polygon.  Default: true
 256///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
 257/// Example(2D):
 258///   path = [
 259///       [-100,100], [0,-50], [100,100], [100,-100], [0,50], [-100,-100]
 260///   ];
 261///   isects = _path_self_intersections(path, closed=true);
 262///   // isects == [[[-33.3333, 0], 0, 0.666667, 4, 0.333333], [[33.3333, 0], 1, 0.333333, 3, 0.666667]]
 263///   stroke(path, closed=true, width=1);
 264///   for (isect=isects) translate(isect[0]) color("blue") sphere(d=10);
 265function _path_self_intersections(path, closed=true, eps=EPSILON) =
 266    let(
 267        path = closed ? list_wrap(path,eps=eps) : path,
 268        plen = len(path)
 269    )
 270    [ for (i = [0:1:plen-3]) let(
 271          a1 = path[i],
 272          a2 = path[i+1], 
 273          seg_normal = unit([-(a2-a1).y, (a2-a1).x],[0,0]),
 274          vals = path*seg_normal,
 275          ref  = a1*seg_normal,
 276            // The value of vals[j]-ref is positive if vertex j is one one side of the
 277            // line [a1,a2] and negative on the other side. Only a segment with opposite
 278            // signs at its two vertices can have an intersection with segment
 279            // [a1,a2]. The variable signals is zero when abs(vals[j]-ref) is less than
 280            // eps and the sign of vals[j]-ref otherwise.  
 281          signals = [for(j=[i+2:1:plen-(i==0 && closed? 2: 1)]) 
 282                        abs(vals[j]-ref) <  eps ? 0 : sign(vals[j]-ref) ]
 283        )
 284        if(max(signals)>=0 && min(signals)<=0 ) // some remaining edge intersects line [a1,a2]
 285        for(j=[i+2:1:plen-(i==0 && closed? 3: 2)])
 286            if( signals[j-i-2]*signals[j-i-1]<=0 ) let( // segm [b1,b2] intersects line [a1,a2]
 287                b1 = path[j],
 288                b2 = path[j+1],
 289                isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) 
 290            )
 291            if (isect 
 292                && isect[1]>=-eps
 293                && isect[1]<= 1+eps
 294                && isect[2]>= -eps 
 295                && isect[2]<= 1+eps)
 296                [isect[0], i, isect[1], j, isect[2]]
 297    ];
 298
 299// Section: Resampling - changing the number of points in a path
 300
 301
 302// Input `data` is a list that sums to an integer. 
 303// Returns rounded version of input data so that every 
 304// entry is rounded to an integer and the sum is the same as
 305// that of the input.  Works by rounding an entry in the list
 306// and passing the rounding error forward to the next entry.
 307// This will generally distribute the error in a uniform manner. 
 308function _sum_preserving_round(data, index=0) =
 309    index == len(data)-1 ? list_set(data, len(data)-1, round(data[len(data)-1])) :
 310    let(
 311        newval = round(data[index]),
 312        error = newval - data[index]
 313    ) _sum_preserving_round(
 314        list_set(data, [index,index+1], [newval, data[index+1]-error]),
 315        index+1
 316    );
 317
 318
 319// Function: subdivide_path()
 320// Synopsis: Subdivides a path to produce a more finely sampled path.
 321// SynTags: Path
 322// Topics: Paths, Path Subdivision
 323// See Also: subdivide_and_slice(), resample_path(), jittered_poly()
 324// Usage:
 325//   newpath = subdivide_path(path, n|refine=|maxlen=, [method=], [closed=], [exact=]);
 326// Description:
 327//   Takes a path as input (closed or open) and subdivides the path to produce a more
 328//   finely sampled path.  You control the subdivision process by using the `maxlen` arg
 329//   to specify a maximum segment length, or by specifying `n` or `refine`, which request
 330//   a certain point count in the output.
 331//   .
 332//   You can specify the point count using the `n` option, where
 333//   you give the number of points you want in the output, or you can use
 334//   the `refine` option, where you specify a resampling factor.  If `refine=3` then
 335//   the number of points would increase by a factor of three, so a four point square would
 336//   have 12 points after subdivision.  With point-count subdivision, the new points can be distributed
 337//   proportional to length (`method="length"`), which is the default, or they can be divided up evenly among all the path segments
 338//   (`method="segment"`).  If the extra points don't fit evenly on the path then the
 339//   algorithm attempts to distribute them as uniformly as possible, but the result may be uneven.
 340//   The `exact` option, which is true by default, requires that the final point count is
 341//   exactly as requested.  For example, if you subdivide a four point square and request `n=13` then one edge will have
 342//   an extra point compared to the others.  
 343//   If you set `exact=false` then the
 344//   algorithm will favor uniformity and the output path may have a different number of
 345//   points than you requested, but the sampling will be uniform.   In our example of the
 346//   square with `n=13`, you will get only 12 points output, with the same number of points on each edge.
 347//   .
 348//   The points are always distributed uniformly on each segment.  The `method="length"` option does
 349//   means that the number of points on a segment is based on its length, but the points are still
 350//   distributed uniformly on each segment, independent of the other segments.  
 351//   With the `"segment"` method you can also give `n` as a vector of counts.  This 
 352//   specifies the desired point count on each segment: with vector valued `n` the `subdivide_path`
 353//   function places `n[i]-1` points on segment `i`.  The reason for the -1 is to avoid
 354//   double counting the endpoints, which are shared by pairs of segments, so that for
 355//   a closed polygon the total number of points will be sum(n).  Note that with an open
 356//   path there is an extra point at the end, so the number of points will be sum(n)+1.
 357//   .
 358//   If you use the `maxlen` option then you specify the maximum length segment allowed in the output.
 359//   Each segment is subdivided into the largest number of segments meeting your requirement.  As above,
 360//   the sampling is uniform on each segment, independent of the other segments.  With the `maxlen` option
 361//   you cannot specify `method` or `exact`.    
 362// Arguments:
 363//   path = path in any dimension or a 1-region
 364//   n = scalar total number of points desired or with `method="segment"` can be a vector requesting `n[i]-1` new points added to segment i.
 365//   ---
 366//   refine = increase total number of points by this factor (Specify only one of n, refine and maxlen)
 367//   maxlen = maximum length segment in the output (Specify only one of n, refine and maxlen)
 368//   closed = set to false if the path is open.  Default: True
 369//   exact = if true return exactly the requested number of points, possibly sacrificing uniformity.  If false, return uniform point sample that may not match the number of points requested.  (Not allowed with maxlen.) Default: true
 370//   method = One of `"length"` or `"segment"`.  If `"length"`, adds vertices in proportion to segment length, so short segments get fewer points.  If `"segment"`, add points evenly among the segments, so all segments get the same number of points.  (Not allowed with maxlen.) Default: `"length"`
 371// Example(2D):
 372//   mypath = subdivide_path(square([2,2],center=true), 12);
 373//   move_copies(mypath)circle(r=.1,$fn=32);
 374// Example(2D):
 375//   mypath = subdivide_path(square([8,2],center=true), 12);
 376//   move_copies(mypath)circle(r=.2,$fn=32);
 377// Example(2D):
 378//   mypath = subdivide_path(square([8,2],center=true), 12, method="segment");
 379//   move_copies(mypath)circle(r=.2,$fn=32);
 380// Example(2D):
 381//   mypath = subdivide_path(square([2,2],center=true), 17, closed=false);
 382//   move_copies(mypath)circle(r=.1,$fn=32);
 383// Example(2D): Specifying different numbers of points on each segment
 384//   mypath = subdivide_path(hexagon(side=2), [2,3,4,5,6,7], method="segment");
 385//   move_copies(mypath)circle(r=.1,$fn=32);
 386// Example(2D): Requested point total is 14 but 15 points output due to extra end point
 387//   mypath = subdivide_path(pentagon(side=2), [3,4,3,4], method="segment", closed=false);
 388//   move_copies(mypath)circle(r=.1,$fn=32);
 389// Example(2D): Since 17 is not divisible by 5, a completely uniform distribution is not possible. 
 390//   mypath = subdivide_path(pentagon(side=2), 17);
 391//   move_copies(mypath)circle(r=.1,$fn=32);
 392// Example(2D): With `exact=false` a uniform distribution, but only 15 points
 393//   mypath = subdivide_path(pentagon(side=2), 17, exact=false);
 394//   move_copies(mypath)circle(r=.1,$fn=32);
 395// Example(2D): With `exact=false` you can also get extra points, here 20 instead of requested 18
 396//   mypath = subdivide_path(pentagon(side=2), 18, exact=false);
 397//   move_copies(mypath)circle(r=.1,$fn=32);
 398// Example(2D): Using refine in this example multiplies the point count by 3 by adding 2 points to each edge
 399//   mypath = subdivide_path(pentagon(side=2), refine=3);
 400//   move_copies(mypath)circle(r=.1,$fn=32);
 401// Example(2D): But note that refine doesn't distribute evenly by segment unless you change the method.  with the default method set to `"length"`, the points are distributed with more on the long segments in this example using refine.  
 402//   mypath = subdivide_path(square([8,2],center=true), refine=3);
 403//   move_copies(mypath)circle(r=.2,$fn=32);
 404// Example(2D): In this example with maxlen, every side gets a different number of new points
 405//   path = [[0,0],[0,4],[10,6],[10,0]];
 406//   spath = subdivide_path(path, maxlen=2, closed=true);
 407//   move_copies(spath) circle(r=.25,$fn=12);
 408// Example(FlatSpin,VPD=15,VPT=[0,0,1.5]): Three-dimensional paths also work
 409//   mypath = subdivide_path([[0,0,0],[2,0,1],[2,3,2]], 12);
 410//   move_copies(mypath)sphere(r=.1,$fn=32);
 411function subdivide_path(path, n, refine, maxlen, closed=true, exact, method) =
 412    let(path = force_path(path))
 413    assert(is_path(path))
 414    assert(num_defined([n,refine,maxlen]),"Must give exactly one of n, refine, and maxlen")
 415    refine==1 || n==len(path) ? path :
 416    is_def(maxlen) ?
 417        assert(is_undef(method), "Cannot give method with maxlen")
 418        assert(is_undef(exact), "Cannot give exact with maxlen")
 419        [
 420         for (p=pair(path,closed))
 421           let(steps = ceil(norm(p[1]-p[0])/maxlen))
 422           each lerpn(p[0], p[1], steps, false),
 423         if (!closed) last(path)
 424        ]               
 425    :
 426    let(
 427        exact = default(exact, true),
 428        method = default(method, "length")
 429    )
 430    assert(method=="length" || method=="segment")
 431    let(
 432        n = !is_undef(n)? n :
 433            !is_undef(refine)? len(path) * refine :
 434            undef
 435    )
 436    assert((is_num(n) && n>0) || is_vector(n),"Parameter n to subdivide_path must be postive number or vector")
 437    let(
 438        count = len(path) - (closed?0:1), 
 439        add_guess = method=="segment"?
 440                       (
 441                          is_list(n)
 442                          ? assert(len(n)==count,"Vector parameter n to subdivide_path has the wrong length")
 443                            add_scalar(n,-1)
 444                          : repeat((n-len(path)) / count, count)
 445                       )
 446                  : // method=="length"
 447                    assert(is_num(n),"Parameter n to subdivide path must be a number when method=\"length\"")
 448                    let(
 449                        path_lens = path_segment_lengths(path,closed),
 450                        add_density = (n - len(path)) / sum(path_lens)
 451                    )
 452                    path_lens * add_density,
 453        add = exact? _sum_preserving_round(add_guess)
 454                   : [for (val=add_guess) round(val)]
 455    )
 456    [
 457        for (i=[0:1:count-1]) 
 458           each lerpn(path[i],select(path,i+1), 1+add[i],endpoint=false),
 459        if (!closed) last(path)
 460    ];
 461
 462
 463
 464
 465// Function: resample_path()
 466// Synopsis: Returns an equidistant set of points along a path.
 467// SynTags: Path
 468// Topics: Paths
 469// See Also: subdivide_path()
 470// Usage:
 471//   newpath = resample_path(path, n|spacing=, [closed=]);
 472// Description:
 473//   Compute a uniform resampling of the input path.  If you specify `n` then the output path will have n
 474//   points spaced uniformly (by linear interpolation along the input path segments).  The only points of the
 475//   input path that are guaranteed to appear in the output path are the starting and ending points, and any
 476//   points that have an angular deflection of at least the number of degrees given in `keep_corners`.
 477//   If you specify `spacing` then the length you give will be rounded to the nearest spacing that gives
 478//   a uniform sampling of the path and the resulting uniformly sampled path is returned.
 479//   Note that because this function operates on a discrete input path the quality of the output depends on
 480//   the sampling of the input.  If you want very accurate output, use a lot of points for the input.
 481// Arguments:
 482//   path = path in any dimension or a 1-region
 483//   n = Number of points in output
 484//   ---
 485//   spacing = Approximate spacing desired
 486//   keep_corners = If given a scalar, path vertices with deflection angle greater than this are preserved in the output.
 487//   closed = set to true if path is closed.  Default: true
 488// Example(2D):  Subsampling lots of points from a smooth curve
 489//   path = xscale(2,circle($fn=250, r=10));
 490//   sampled = resample_path(path, 16);
 491//   stroke(path);
 492//   color("red")move_copies(sampled) circle($fn=16);
 493// Example(2D): Specified spacing is rounded to make a uniform sampling
 494//   path = xscale(2,circle($fn=250, r=10));
 495//   sampled = resample_path(path, spacing=17);
 496//   stroke(path);
 497//   color("red")move_copies(sampled) circle($fn=16);
 498// Example(2D): Notice that the corners are excluded.
 499//   path = square(20);
 500//   sampled = resample_path(path, spacing=6);
 501//   stroke(path,closed=true);
 502//   color("red")move_copies(sampled) circle($fn=16);
 503// Example(2D): Forcing preservation of corners.
 504//   path = square(20);
 505//   sampled = resample_path(path, spacing=6, keep_corners=90);
 506//   stroke(path,closed=true);
 507//   color("red")move_copies(sampled) circle($fn=16);
 508// Example(2D): Closed set to false
 509//   path = square(20);
 510//   sampled = resample_path(path, spacing=6,closed=false);
 511//   stroke(path);
 512//   color("red")move_copies(sampled) circle($fn=16);
 513
 514function resample_path(path, n, spacing, keep_corners, closed=true) =
 515    let(path = force_path(path))
 516    assert(is_path(path))
 517    assert(num_defined([n,spacing])==1,"Must define exactly one of n and spacing")
 518    assert(n==undef || (is_integer(n) && n>0))
 519    assert(spacing==undef || (is_finite(spacing) && spacing>0))
 520    assert(is_bool(closed))
 521    let(
 522        corners = is_undef(keep_corners)
 523          ? [0, len(path)-(closed?0:1)]
 524          : [
 525                0,
 526                for (i = [1:1:len(path)-(closed?1:2)])
 527                    let( ang = abs(modang(vector_angle(select(path,i-1,i+1))-180)) )
 528                    if (ang >= keep_corners) i,
 529                len(path)-(closed?0:1),
 530            ],
 531        pcnt = len(path),
 532        plen = path_length(path, closed=closed),
 533        subpaths = [ for (p = pair(corners)) [for(i = [p.x:1:p.y]) path[i%pcnt]] ],
 534        n = is_undef(n)? undef : closed? n+1 : n
 535    )
 536    assert(n==undef || n >= len(corners), "There are nore than `n=` corners whose angle is greater than `keep_corners=`.")
 537    let(
 538        lens = [for (subpath = subpaths) path_length(subpath)],
 539        part_ns = is_undef(n)
 540          ? [for (i=idx(subpaths)) max(1,round(lens[i]/spacing)-1)]
 541          : let(
 542                ccnt = len(corners),
 543                parts = [for (l=lens) (n-ccnt) * l/plen]
 544            )
 545            _sum_preserving_round(parts),
 546        out = [
 547            for (i = idx(subpaths))
 548                let(
 549                    subpath = subpaths[i],
 550                    splen = lens[i],
 551                    pn = part_ns[i] + 1,
 552                    distlist = lerpn(0, splen, pn, false),
 553                    cuts = path_cut_points(subpath, distlist, closed=false)
 554                )
 555                each column(cuts,0),
 556            if (!closed) last(path)
 557        ]
 558    ) out;
 559
 560
 561// Section: Path Geometry
 562
 563// Function: is_path_simple()
 564// Synopsis: Returns true if a path has no self intersections.
 565// Topics: Paths
 566// See Also: is_path()
 567// Usage:
 568//   bool = is_path_simple(path, [closed], [eps]);
 569// Description:
 570//   Returns true if the given 2D path is simple, meaning that it has no self-intersections.
 571//   Repeated points are not considered self-intersections: a path with such points can
 572//   still be simple.  
 573//   If closed is set to true then treat the path as a polygon.
 574// Arguments:
 575//   path = 2D path or 1-region
 576//   closed = set to true to treat path as a polygon.  Default: false
 577//   eps = Epsilon error value used for determine if points coincide.  Default: `EPSILON` (1e-9)
 578function is_path_simple(path, closed, eps=EPSILON) =
 579    is_1region(path) ? is_path_simple(path[0], default(closed,true), eps) :
 580    let(closed=default(closed,false))
 581    assert(is_path(path, 2),"Must give a 2D path")
 582    assert(is_bool(closed))
 583    let(
 584        path = deduplicate(path,closed=closed,eps=eps)
 585    )
 586    // check for path reversals
 587    [for(i=[0:1:len(path)-(closed?2:3)])
 588         let(v1=path[i+1]-path[i],
 589             v2=select(path,i+2)-path[i+1],
 590             normv1 = norm(v1),
 591             normv2 = norm(v2)
 592             )
 593         if (approx(v1*v2/normv1/normv2,-1)) 1
 594    ]  == [] 
 595    &&
 596    _path_self_intersections(path,closed=closed,eps=eps) == [];
 597
 598
 599// Function: path_closest_point()
 600// Synopsis: Returns the closest place on a path to a given point.
 601// Topics: Paths
 602// See Also: point_line_distance(), line_closest_point()
 603// Usage:
 604//   index_pt = path_closest_point(path, pt);
 605// Description:
 606//   Finds the closest path segment, and point on that segment to the given point.
 607//   Returns `[SEGNUM, POINT]`
 608// Arguments:
 609//   path = Path of any dimension or a 1-region.
 610//   pt = The point to find the closest point to.
 611//   closed = If true, the path is considered closed.
 612// Example(2D):
 613//   path = circle(d=100,$fn=6);
 614//   pt = [20,10];
 615//   closest = path_closest_point(path, pt);
 616//   stroke(path, closed=true);
 617//   color("blue") translate(pt) circle(d=3, $fn=12);
 618//   color("red") translate(closest[1]) circle(d=3, $fn=12);
 619function path_closest_point(path, pt, closed=true) =
 620    let(path = force_path(path))
 621    assert(is_path(path), "Input must be a path")
 622    assert(is_vector(pt, len(path[0])), "Input pt must be a compatible vector")
 623    assert(is_bool(closed))
 624    let(
 625        pts = [for (seg=pair(path,closed)) line_closest_point(seg,pt,SEGMENT)],
 626        dists = [for (p=pts) norm(p-pt)],
 627        min_seg = min_index(dists)
 628    ) [min_seg, pts[min_seg]];
 629
 630
 631// Function: path_tangents()
 632// Synopsis: Returns tangent vectors for each point along a path.
 633// Topics: Paths
 634// See Also: path_normals()
 635// Usage:
 636//   tangs = path_tangents(path, [closed], [uniform]);
 637// Description:
 638//   Compute the tangent vector to the input path.  The derivative approximation is described in deriv().
 639//   The returns vectors will be normalized to length 1.  If any derivatives are zero then
 640//   the function fails with an error.  If you set `uniform` to false then the sampling is
 641//   assumed to be non-uniform and the derivative is computed with adjustments to produce corrected
 642//   values.
 643// Arguments:
 644//   path = path of any dimension or a 1-region
 645//   closed = set to true of the path is closed.  Default: false
 646//   uniform = set to false to correct for non-uniform sampling.  Default: true
 647// Example(2D): A shape with non-uniform sampling gives distorted derivatives that may be undesirable.  Note that derivatives tilt towards the long edges of the rectangle.  
 648//   rect = square([10,3]);
 649//   tangents = path_tangents(rect,closed=true);
 650//   stroke(rect,closed=true, width=0.25);
 651//   color("purple")
 652//       for(i=[0:len(tangents)-1])
 653//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.25, endcap2="arrow2");
 654// Example(2D): Setting uniform to false corrects the distorted derivatives for this example:
 655//   rect = square([10,3]);
 656//   tangents = path_tangents(rect,closed=true,uniform=false);
 657//   stroke(rect,closed=true, width=0.25);
 658//   color("purple")
 659//       for(i=[0:len(tangents)-1])
 660//           stroke([rect[i]-tangents[i], rect[i]+tangents[i]],width=.25, endcap2="arrow2");
 661function path_tangents(path, closed, uniform=true) =
 662    is_1region(path) ? path_tangents(path[0], default(closed,true), uniform) :
 663    let(closed=default(closed,false))
 664    assert(is_bool(closed))
 665    assert(is_path(path))
 666    !uniform ? [for(t=deriv(path,closed=closed, h=path_segment_lengths(path,closed))) unit(t)]
 667             : [for(t=deriv(path,closed=closed)) unit(t)];
 668
 669
 670// Function: path_normals()
 671// Synopsis: Returns normal vectors for each point along a path.
 672// Topics: Paths
 673// See Also: path_tangents()
 674// Usage:
 675//   norms = path_normals(path, [tangents], [closed]);
 676// Description:
 677//   Compute the normal vector to the input path.  This vector is perpendicular to the
 678//   path tangent and lies in the plane of the curve.  For 3d paths we define the plane of the curve
 679//   at path point i to be the plane defined by point i and its two neighbors.  At the endpoints of open paths
 680//   we use the three end points.  For 3d paths the computed normal is the one lying in this plane that points
 681//   towards the center of curvature at that path point.  For 2d paths, which lie in the xy plane, the normal
 682//   is the path pointing to the right of the direction the path is traveling.  If points are collinear then
 683//   a 3d path has no center of curvature, and hence the 
 684//   normal is not uniquely defined.  In this case the function issues an error.
 685//   For 2d paths the plane is always defined so the normal fails to exist only
 686//   when the derivative is zero (in the case of repeated points).
 687// Arguments:
 688//   path = 2D or 3D path or a 1-region
 689//   tangents = path tangents optionally supplied
 690//   closed = if true path is treated as a polygon.  Default: false
 691function path_normals(path, tangents, closed) =
 692    is_1region(path) ? path_normals(path[0], tangents, default(closed,true)) :
 693    let(closed=default(closed,false))
 694    assert(is_path(path,[2,3]))
 695    assert(is_bool(closed))
 696    let(
 697         tangents = default(tangents, path_tangents(path,closed)),
 698         dim=len(path[0])
 699    )
 700    assert(is_path(tangents) && len(tangents[0])==dim,"Dimensions of path and tangents must match")
 701    [
 702     for(i=idx(path))
 703         let(
 704             pts = i==0 ? (closed? select(path,-1,1) : select(path,0,2))
 705                 : i==len(path)-1 ? (closed? select(path,i-1,i+1) : select(path,i-2,i))
 706                 : select(path,i-1,i+1)
 707        )
 708        dim == 2 ? [tangents[i].y,-tangents[i].x]
 709                 : let( v=cross(cross(pts[1]-pts[0], pts[2]-pts[0]),tangents[i]))
 710                   assert(norm(v)>EPSILON, "3D path contains collinear points")
 711                   unit(v)
 712    ];
 713
 714
 715// Function: path_curvature()
 716// Synopsis: Returns the estimated numerical curvature of the path.
 717// Topics: Paths
 718// See Also: path_tangents(), path_normals(), path_torsion()
 719// Usage:
 720//   curvs = path_curvature(path, [closed]);
 721// Description:
 722//   Numerically estimate the curvature of the path (in any dimension).
 723// Arguments:
 724//   path = path in any dimension or a 1-region
 725//   closed = if true then treat the path as a polygon.  Default: false
 726function path_curvature(path, closed) =
 727    is_1region(path) ? path_curvature(path[0], default(closed,true)) :
 728    let(closed=default(closed,false))
 729    assert(is_bool(closed))
 730    assert(is_path(path))
 731    let( 
 732        d1 = deriv(path, closed=closed),
 733        d2 = deriv2(path, closed=closed)
 734    ) [
 735        for(i=idx(path))
 736        sqrt(
 737            sqr(norm(d1[i])*norm(d2[i])) -
 738            sqr(d1[i]*d2[i])
 739        ) / pow(norm(d1[i]),3)
 740    ];
 741
 742
 743// Function: path_torsion()
 744// Synopsis: Returns the estimated numerical torsion of the path.
 745// Topics: Paths
 746// See Also: path_tangents(), path_normals(), path_curvature()
 747// Usage:
 748//   torsions = path_torsion(path, [closed]);
 749// Description:
 750//   Numerically estimate the torsion of a 3d path.
 751// Arguments:
 752//   path = 3D path
 753//   closed = if true then treat path as a polygon.  Default: false
 754function path_torsion(path, closed=false) =
 755    assert(is_path(path,3), "Input path must be a 3d path")
 756    assert(is_bool(closed))
 757    let(
 758        d1 = deriv(path,closed=closed),
 759        d2 = deriv2(path,closed=closed),
 760        d3 = deriv3(path,closed=closed)
 761    ) [
 762        for (i=idx(path)) let(
 763            crossterm = cross(d1[i],d2[i])
 764        ) crossterm * d3[i] / sqr(norm(crossterm))
 765    ];
 766
 767
 768// Section: Breaking paths up into subpaths
 769
 770
 771
 772// Function: path_cut()
 773// Synopsis: Cuts a path into subpaths at various points.
 774// SynTags: PathList
 775// Topics: Paths, Path Subdivision
 776// See Also: split_path_at_self_crossings(), path_cut_points()
 777// Usage:
 778//   path_list = path_cut(path, cutdist, [closed]);
 779// Description:
 780//   Given a list of distances in `cutdist`, cut the path into
 781//   subpaths at those lengths, returning a list of paths.
 782//   If the input path is closed then the final path will include the
 783//   original starting point.  The list of cut distances must be
 784//   in ascending order and should not include the endpoints: 0 
 785//   or len(path).  If you repeat a distance you will get an
 786//   empty list in that position in the output.  If you give an
 787//   empty cutdist array you will get the input path as output
 788//   (without the final vertex doubled in the case of a closed path).
 789// Arguments:
 790//   path = path of any dimension or a 1-region
 791//   cutdist = Distance or list of distances where path is cut
 792//   closed = If true, treat the path as a closed polygon.  Default: false
 793// Example(2D,NoAxes):
 794//   path = circle(d=100);
 795//   segs = path_cut(path, [50, 200], closed=true);
 796//   rainbow(segs) stroke($item, endcaps="butt", width=3);
 797function path_cut(path,cutdist,closed) =
 798  is_num(cutdist) ? path_cut(path,[cutdist],closed) :
 799  is_1region(path) ? path_cut(path[0], cutdist, default(closed,true)):
 800  let(closed=default(closed,false))
 801  assert(is_bool(closed))
 802  assert(is_vector(cutdist))
 803  assert(last(cutdist)<path_length(path,closed=closed)-EPSILON,"Cut distances must be smaller than the path length")
 804  assert(cutdist[0]>EPSILON, "Cut distances must be strictly positive")
 805  let(
 806      cutlist = path_cut_points(path,cutdist,closed=closed)
 807  )
 808  _path_cut_getpaths(path, cutlist, closed);
 809
 810
 811function _path_cut_getpaths(path, cutlist, closed) =
 812  let(
 813      cuts = len(cutlist)
 814  )
 815  [
 816      [ each list_head(path,cutlist[0][1]-1),
 817        if (!approx(cutlist[0][0], path[cutlist[0][1]-1])) cutlist[0][0]
 818      ],
 819      for(i=[0:1:cuts-2])
 820          cutlist[i][0]==cutlist[i+1][0] && cutlist[i][1]==cutlist[i+1][1] ? []
 821          :
 822          [ if (!approx(cutlist[i][0], select(path,cutlist[i][1]))) cutlist[i][0],
 823            each slice(path, cutlist[i][1], cutlist[i+1][1]-1),
 824            if (!approx(cutlist[i+1][0], select(path,cutlist[i+1][1]-1))) cutlist[i+1][0],
 825          ],
 826      [
 827        if (!approx(cutlist[cuts-1][0], select(path,cutlist[cuts-1][1]))) cutlist[cuts-1][0],
 828        each select(path,cutlist[cuts-1][1],closed ? 0 : -1)
 829      ]
 830  ];
 831
 832
 833
 834// Function: path_cut_points()
 835// Synopsis: Returns a list of cut points at a list of distances from the first point in a path.
 836// Topics: Paths, Path Subdivision
 837// See Also: path_cut(), split_path_at_self_crossings()
 838// Usage:
 839//   cuts = path_cut_points(path, cutdist, [closed=], [direction=]);
 840//
 841// Description:
 842//   Cuts a path at a list of distances from the first point in the path.  Returns a list of the cut
 843//   points and indices of the next point in the path after that point.  So for example, a return
 844//   value entry of [[2,3], 5] means that the cut point was [2,3] and the next point on the path after
 845//   this point is path[5].  If the path is too short then path_cut_points returns undef.  If you set
 846//   `direction` to true then `path_cut_points` will also return the tangent vector to the path and a normal
 847//   vector to the path.  It tries to find a normal vector that is coplanar to the path near the cut
 848//   point.  If this fails it will return a normal vector parallel to the xy plane.  The output with
 849//   direction vectors will be `[point, next_index, tangent, normal]`.
 850//   .
 851//   If you give the very last point of the path as a cut point then the returned index will be
 852//   one larger than the last index (so it will not be a valid index).  If you use the closed
 853//   option then the returned index will be equal to the path length for cuts along the closing
 854//   path segment, and if you give a point equal to the path length you will get an
 855//   index of len(path)+1 for the index.  
 856//
 857// Arguments:
 858//   path = path to cut
 859//   cutdist = distances where the path should be cut (a list) or a scalar single distance
 860//   ---
 861//   closed = set to true if the curve is closed.  Default: false
 862//   direction = set to true to return direction vectors.  Default: false
 863//
 864// Example(NORENDER):
 865//   square=[[0,0],[1,0],[1,1],[0,1]];
 866//   path_cut_points(square, [.5,1.5,2.5]);   // Returns [[[0.5, 0], 1], [[1, 0.5], 2], [[0.5, 1], 3]]
 867//   path_cut_points(square, [0,1,2,3]);      // Returns [[[0, 0], 1], [[1, 0], 2], [[1, 1], 3], [[0, 1], 4]]
 868//   path_cut_points(square, [0,0.8,1.6,2.4,3.2], closed=true);  // Returns [[[0, 0], 1], [[0.8, 0], 1], [[1, 0.6], 2], [[0.6, 1], 3], [[0, 0.8], 4]]
 869//   path_cut_points(square, [0,0.8,1.6,2.4,3.2]);               // Returns [[[0, 0], 1], [[0.8, 0], 1], [[1, 0.6], 2], [[0.6, 1], 3], undef]
 870function path_cut_points(path, cutdist, closed=false, direction=false) =
 871    let(long_enough = len(path) >= (closed ? 3 : 2))
 872    assert(long_enough,len(path)<2 ? "Two points needed to define a path" : "Closed path must include three points")
 873    is_num(cutdist) ? path_cut_points(path, [cutdist],closed, direction)[0] :
 874    assert(is_vector(cutdist))
 875    assert(is_increasing(cutdist), "Cut distances must be an increasing list")
 876    let(cuts = path_cut_points_recurse(path,cutdist,closed))
 877    !direction
 878       ? cuts
 879       : let(
 880             dir = _path_cuts_dir(path, cuts, closed),
 881             normals = _path_cuts_normals(path, cuts, dir, closed)
 882         )
 883         hstack(cuts, list_to_matrix(dir,1), list_to_matrix(normals,1));
 884
 885// Main recursive path cut function
 886function path_cut_points_recurse(path, dists, closed=false, pind=0, dtotal=0, dind=0, result=[]) =
 887    dind == len(dists) ? result :
 888    let(
 889        lastpt = len(result)==0? [] : last(result)[0],       // location of last cut point
 890        dpartial = len(result)==0? 0 : norm(lastpt-select(path,pind)),  // remaining length in segment
 891        nextpoint = dists[dind] < dpartial+dtotal  // Do we have enough length left on the current segment?
 892           ? [lerp(lastpt,select(path,pind),(dists[dind]-dtotal)/dpartial),pind] 
 893           : _path_cut_single(path, dists[dind]-dtotal-dpartial, closed, pind)
 894    ) 
 895    path_cut_points_recurse(path, dists, closed, nextpoint[1], dists[dind],dind+1, concat(result, [nextpoint]));
 896
 897
 898// Search for a single cut point in the path
 899function _path_cut_single(path, dist, closed=false, ind=0, eps=1e-7) =
 900    // If we get to the very end of the path (ind is last point or wraparound for closed case) then
 901    // check if we are within epsilon of the final path point.  If not we're out of path, so we fail
 902    ind==len(path)-(closed?0:1) ?
 903       assert(dist<eps,"Path is too short for specified cut distance")
 904       [select(path,ind),ind+1]
 905    :let(d = norm(path[ind]-select(path,ind+1))) d > dist ?
 906        [lerp(path[ind],select(path,ind+1),dist/d), ind+1] :
 907        _path_cut_single(path, dist-d,closed, ind+1, eps);
 908
 909// Find normal directions to the path, coplanar to local part of the path
 910// Or return a vector parallel to the x-y plane if the above fails
 911function _path_cuts_normals(path, cuts, dirs, closed=false) =
 912    [for(i=[0:len(cuts)-1])
 913        len(path[0])==2? [-dirs[i].y, dirs[i].x]
 914          : 
 915            let(
 916                plane = len(path)<3 ? undef :
 917                let(start = max(min(cuts[i][1],len(path)-1),2)) _path_plane(path, start, start-2)
 918            )
 919            plane==undef?
 920                ( dirs[i].x==0 && dirs[i].y==0 ? [1,0,0]  // If it's z direction return x vector
 921                                               : unit([-dirs[i].y, dirs[i].x,0])) // otherwise perpendicular to projection
 922                : unit(cross(dirs[i],cross(plane[0],plane[1])))
 923    ];
 924
 925// Scan from the specified point (ind) to find a noncoplanar triple to use
 926// to define the plane of the path.
 927function _path_plane(path, ind, i,closed) =
 928    i<(closed?-1:0) ? undef :
 929    !is_collinear(path[ind],path[ind-1], select(path,i))?
 930        [select(path,i)-path[ind-1],path[ind]-path[ind-1]] :
 931        _path_plane(path, ind, i-1);
 932
 933// Find the direction of the path at the cut points
 934function _path_cuts_dir(path, cuts, closed=false, eps=1e-2) =
 935    [for(ind=[0:len(cuts)-1])
 936        let(
 937            zeros = path[0]*0,
 938            nextind = cuts[ind][1],
 939            nextpath = unit(select(path, nextind+1)-select(path, nextind),zeros),
 940            thispath = unit(select(path, nextind) - select(path,nextind-1),zeros),
 941            lastpath = unit(select(path,nextind-1) - select(path, nextind-2),zeros),
 942            nextdir =
 943                nextind==len(path) && !closed? lastpath :
 944                (nextind<=len(path)-2 || closed) && approx(cuts[ind][0], path[nextind],eps)
 945                   ? unit(nextpath+thispath)
 946              : (nextind>1 || closed) && approx(cuts[ind][0],select(path,nextind-1),eps)
 947                   ? unit(thispath+lastpath)
 948              :  thispath
 949        ) nextdir
 950    ];
 951
 952
 953// internal function
 954// converts pathcut output form to a [segment, u]
 955// form list that works withi path_select
 956function _cut_to_seg_u_form(pathcut, path, closed) =
 957  let(lastind = len(path) - (closed?0:1))
 958  [for(entry=pathcut)
 959    entry[1] > lastind ? [lastind,0] :
 960    let(
 961        a = path[entry[1]-1],
 962        b = path[entry[1]],
 963        c = entry[0],
 964        i = max_index(v_abs(b-a)),
 965        factor = (c[i]-a[i])/(b[i]-a[i])
 966    )
 967    [entry[1]-1,factor]
 968  ];
 969
 970
 971
 972// Function: split_path_at_self_crossings()
 973// Synopsis: Split a 2D path wherever it crosses itself.
 974// SynTags: PathList
 975// Topics: Paths, Path Subdivision
 976// See Also: path_cut(), path_cut_points()
 977// Usage:
 978//   paths = split_path_at_self_crossings(path, [closed], [eps]);
 979// Description:
 980//   Splits a 2D path into sub-paths wherever the original path crosses itself.
 981//   Splits may occur mid-segment, so new vertices will be created at the intersection points.
 982//   Returns a list of the resulting subpaths.  
 983// Arguments:
 984//   path = A 2D path or a 1-region.
 985//   closed = If true, treat path as a closed polygon.  Default: true
 986//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 987// Example(2D,NoAxes):
 988//   path = [ [-100,100], [0,-50], [100,100], [100,-100], [0,50], [-100,-100] ];
 989//   paths = split_path_at_self_crossings(path);
 990//   rainbow(paths) stroke($item, closed=false, width=3);
 991function split_path_at_self_crossings(path, closed=true, eps=EPSILON) =
 992    let(path = force_path(path))
 993    assert(is_path(path,2), "Must give a 2D path")
 994    assert(is_bool(closed))
 995    let(
 996        path = list_unwrap(path, eps=eps),
 997        isects = deduplicate(
 998            eps=eps,
 999            concat(
1000                [[0, 0]],
1001                sort([
1002                    for (
1003                        a = _path_self_intersections(path, closed=closed, eps=eps),
1004                        ss = [ [a[1],a[2]], [a[3],a[4]] ]
1005                    ) if (ss[0] != undef) ss
1006                ]),
1007                [[len(path)-(closed?1:2), 1]]
1008            )
1009        )
1010    ) [
1011        for (p = pair(isects))
1012            let(
1013                s1 = p[0][0],
1014                u1 = p[0][1],
1015                s2 = p[1][0],
1016                u2 = p[1][1],
1017                section = _path_select(path, s1, u1, s2, u2, closed=closed),
1018                outpath = deduplicate(eps=eps, section)
1019            )
1020            if (len(outpath)>1) outpath
1021    ];
1022
1023
1024function _tag_self_crossing_subpaths(path, nonzero, closed=true, eps=EPSILON) =
1025    let(
1026        subpaths = split_path_at_self_crossings(
1027            path, closed=true, eps=eps
1028        )
1029    ) [
1030        for (subpath = subpaths) let(
1031            seg = select(subpath,0,1),
1032            mp = mean(seg),
1033            n = line_normal(seg) / 2048,
1034            p1 = mp + n,
1035            p2 = mp - n,
1036            p1in = point_in_polygon(p1, path, nonzero=nonzero) >= 0,
1037            p2in = point_in_polygon(p2, path, nonzero=nonzero) >= 0,
1038            tag = (p1in && p2in)? "I" : "O"
1039        ) [tag, subpath]
1040    ];
1041
1042
1043// Function: polygon_parts()
1044// Synopsis: Parses a self-intersecting polygon into a list of non-intersecting polygons.
1045// SynTags: PathList
1046// Topics: Paths, Polygons
1047// See Also: split_path_at_self_crossings(), path_cut(), path_cut_points()
1048// Usage:
1049//   splitpolys = polygon_parts(poly, [nonzero], [eps]);
1050// Description:
1051//   Given a possibly self-intersecting 2d polygon, constructs a representation of the original polygon as a list of
1052//   non-intersecting simple polygons.  If nonzero is set to true then it uses the nonzero method for defining polygon membership.
1053//   For simple cases, such as the pentagram, this will produce the outer perimeter of a self-intersecting polygon.  
1054// Arguments:
1055//   poly = a 2D polygon or 1-region
1056//   nonzero = If true use the nonzero method for checking if a point is in a polygon.  Otherwise use the even-odd method.  Default: false
1057//   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1058// Example(2D,NoAxes):  This cross-crossing polygon breaks up into its 3 components (regardless of the value of nonzero).
1059//   poly = [
1060//       [-100,100], [0,-50], [100,100],
1061//       [100,-100], [0,50], [-100,-100]
1062//   ];
1063//   splitpolys = polygon_parts(poly);
1064//   rainbow(splitpolys) stroke($item, closed=true, width=3);
1065// Example(2D,NoAxes): With nonzero=false you get even-odd mode which matches OpenSCAD, so the pentagram breaks apart into its five points.
1066//   pentagram = turtle(["move",100,"left",144], repeat=4);
1067//   left(100)polygon(pentagram);
1068//   rainbow(polygon_parts(pentagram,nonzero=false))
1069//     stroke($item,closed=true,width=2.5);
1070// Example(2D,NoAxes): With nonzero=true you get only the outer perimeter.  You can use this to create the polygon using the nonzero method, which is not supported by OpenSCAD.
1071//   pentagram = turtle(["move",100,"left",144], repeat=4);
1072//   outside = polygon_parts(pentagram,nonzero=true);
1073//   left(100)region(outside);
1074//   rainbow(outside)
1075//     stroke($item,closed=true,width=2.5);
1076// Example(2D,NoAxes): 
1077//   N=12;
1078//   ang=360/N;
1079//   sr=10;
1080//   poly = turtle(["angle", 90+ang/2,
1081//                  "move", sr, "left",
1082//                  "move", 2*sr*sin(ang/2), "left",
1083//                  "repeat", 4,
1084//                     ["move", 2*sr, "left",
1085//                      "move", 2*sr*sin(ang/2), "left"],
1086//                  "move", sr]);
1087//   stroke(poly, width=.3);
1088//   right(20)rainbow(polygon_parts(poly)) polygon($item);
1089// Example(2D,NoAxes): overlapping poly segments disappear
1090//   poly = [[0,0], [10,0], [10,10], [0,10],[0,20], [20,10],[10,10], [0,10],[0,0]];
1091//   stroke(poly,width=0.3);
1092//   right(22)stroke(polygon_parts(poly)[0], width=0.3, closed=true);
1093// Example(2D,NoAxes): Poly segments disappear outside as well
1094//   poly = turtle(["repeat", 3, ["move", 17, "left", "move", 10, "left", "move", 7, "left", "move", 10, "left"]]);
1095//   back(2)stroke(poly,width=.5);
1096//   fwd(12)rainbow(polygon_parts(poly)) stroke($item, closed=true, width=0.5);
1097// Example(2D,NoAxes):  This shape has six components
1098//   poly = turtle(["repeat", 3, ["move", 15, "left", "move", 7, "left", "move", 10, "left", "move", 17, "left"]]);
1099//   polygon(poly);
1100//   right(22)rainbow(polygon_parts(poly)) polygon($item);
1101// Example(2D,NoAxes): When the loops of the shape overlap then nonzero gives a different result than the even-odd method.
1102//   poly = turtle(["repeat", 3, ["move", 15, "left", "move", 7, "left", "move", 10, "left", "move", 10, "left"]]);
1103//   polygon(poly);
1104//   right(27)rainbow(polygon_parts(poly)) polygon($item);
1105//   move([16,-14])rainbow(polygon_parts(poly,nonzero=true)) polygon($item);
1106function polygon_parts(poly, nonzero=false, eps=EPSILON) =
1107    let(poly = force_path(poly))
1108    assert(is_path(poly,2), "Must give 2D polygon")
1109    assert(is_bool(nonzero))    
1110    let(
1111        poly = list_unwrap(poly, eps=eps),
1112        tagged = _tag_self_crossing_subpaths(poly, nonzero=nonzero, closed=true, eps=eps),
1113        kept = [for (sub = tagged) if(sub[0] == "O") sub[1]],
1114        outregion = _assemble_path_fragments(kept, eps=eps)
1115    ) outregion;
1116
1117
1118function _extreme_angle_fragment(seg, fragments, rightmost=true, eps=EPSILON) =
1119    !fragments? [undef, []] :
1120    let(
1121        delta = seg[1] - seg[0],
1122        segang = atan2(delta.y,delta.x),
1123        frags = [
1124            for (i = idx(fragments)) let(
1125                fragment = fragments[i],
1126                fwdmatch = approx(seg[1], fragment[0], eps=eps),
1127                bakmatch =  approx(seg[1], last(fragment), eps=eps)
1128            ) [
1129                fwdmatch,
1130                bakmatch,
1131                bakmatch? reverse(fragment) : fragment
1132            ]
1133        ],
1134        angs = [
1135            for (frag = frags)
1136                (frag[0] || frag[1])? let(
1137                    delta2 = frag[2][1] - frag[2][0],
1138                    segang2 = atan2(delta2.y, delta2.x)
1139                ) modang(segang2 - segang) : (
1140                    rightmost? 999 : -999
1141                )
1142        ],
1143        fi = rightmost? min_index(angs) : max_index(angs)
1144    ) abs(angs[fi]) > 360? [undef, fragments] : let(
1145        remainder = [for (i=idx(fragments)) if (i!=fi) fragments[i]],
1146        frag = frags[fi],
1147        foundfrag = frag[2]
1148    ) [foundfrag, remainder];
1149
1150
1151/// Internal Function: _assemble_a_path_from_fragments()
1152/// Usage:
1153///   _assemble_a_path_from_fragments(subpaths);
1154/// Description:
1155///   Given a list of paths, assembles them together into one complete closed polygon path, and
1156///   remainder fragments.  Returns [PATH, FRAGMENTS] where FRAGMENTS is the list of remaining
1157///   unused path fragments.
1158/// Arguments:
1159///   fragments = List of paths to be assembled into complete polygons.
1160///   rightmost = If true, assemble paths using rightmost turns. Leftmost if false.
1161///   startfrag = The fragment to start with.  Default: 0
1162///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1163function _assemble_a_path_from_fragments(fragments, rightmost=true, startfrag=0, eps=EPSILON) =
1164    len(fragments)==0? [[],[]] :
1165    len(fragments)==1? [fragments[0],[]] :
1166    let(
1167        path = fragments[startfrag],
1168        newfrags = [for (i=idx(fragments)) if (i!=startfrag) fragments[i]]
1169    ) are_ends_equal(path, eps=eps)? (
1170        // starting fragment is already closed
1171        [path, newfrags]
1172    ) : let(
1173        // Find rightmost/leftmost continuation fragment
1174        seg = select(path,-2,-1),
1175        extrema = _extreme_angle_fragment(seg=seg, fragments=newfrags, rightmost=rightmost, eps=eps),
1176        foundfrag = extrema[0],
1177        remainder = extrema[1]
1178    ) is_undef(foundfrag)? (
1179        // No remaining fragments connect!  INCOMPLETE PATH!
1180        // Treat it as complete.
1181        [path, remainder]
1182    ) : are_ends_equal(foundfrag, eps=eps)? (
1183        // Found fragment is already closed
1184        [foundfrag, concat([path], remainder)]
1185    ) : let(
1186        fragend = last(foundfrag),
1187        hits = [for (i = idx(path,e=-2)) if(approx(path[i],fragend,eps=eps)) i]
1188    ) hits? (
1189        let(
1190            // Found fragment intersects with initial path
1191            hitidx = last(hits),
1192            newpath = list_head(path,hitidx),
1193            newfrags = concat(len(newpath)>1? [newpath] : [], remainder),
1194            outpath = concat(slice(path,hitidx,-2), foundfrag)
1195        )
1196        [outpath, newfrags]
1197    ) : let(
1198        // Path still incomplete.  Continue building it.
1199        newpath = concat(path, list_tail(foundfrag)),
1200        newfrags = concat([newpath], remainder)
1201    )
1202    _assemble_a_path_from_fragments(
1203        fragments=newfrags,
1204        rightmost=rightmost,
1205        eps=eps
1206    );
1207
1208
1209/// Internal Function: _assemble_path_fragments()
1210/// Usage:
1211///   _assemble_path_fragments(subpaths);
1212/// Description:
1213///   Given a list of paths, assembles them together into complete closed polygon paths if it can.
1214///   Polygons with area < eps will be discarded and not returned.  
1215/// Arguments:
1216///   fragments = List of paths to be assembled into complete polygons.
1217///   eps = The epsilon error value to determine whether two points coincide.  Default: `EPSILON` (1e-9)
1218function _assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
1219    len(fragments)==0? _finished :
1220    let(
1221        minxidx = min_index([
1222            for (frag=fragments) min(column(frag,0))
1223        ]),
1224        result_l = _assemble_a_path_from_fragments(
1225            fragments=fragments,
1226            startfrag=minxidx,
1227            rightmost=false,
1228            eps=eps
1229        ),
1230        result_r = _assemble_a_path_from_fragments(
1231            fragments=fragments,
1232            startfrag=minxidx,
1233            rightmost=true,
1234            eps=eps
1235        ),
1236        l_area = abs(polygon_area(result_l[0])),
1237        r_area = abs(polygon_area(result_r[0])),
1238        result = l_area < r_area? result_l : result_r,
1239        newpath = list_unwrap(result[0]),
1240        remainder = result[1],
1241        finished = min(l_area,r_area)<eps ? _finished : concat(_finished, [newpath])
1242    ) _assemble_path_fragments(
1243        fragments=remainder,
1244        eps=eps,
1245        _finished=finished
1246    );
1247
1248
1249
1250// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap